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Abstract. We extend the framework of the finite volume method to dispersive unidi- 
rectional water wave propagation in one space dimension. In particular we consider a 
KdV-BBM type equation. Explicit and IMEX Runge-Kutta type methods are used for 
time discretizations. The fully discrete schemes are validated by direct comparisons to 
analytic solutions. Invariants conservation properties are also studied. Main applications 
include important nonlinear phenomena such as dispersive shock wave formation, solitary 
waves and their various interactions. 



Water wave modeling is a complicated process and usually leads to models which are 
hard to analyze mathematically as well as to solve numerically. Under certain simplifying 
assumptions approximate models are obtained, e.g. the KdV equation [28], the BBM 
equation [4] and Boussinesq systems [11, 40, 8]. All these models assume the wave to be 
weakly nonlinear and weakly dispersive, propagating mainly in one space direction. These 
approximate models consider mainly unidirectional or bidirectional wave propagation on 
flat or complex bathymetries. 

In this paper we study the application of some finite volume schemes to a scalar non- 
linear dispersive partial differential equation modeling unidirectional wave propagation. 
Specifically, we consider the KdV-BBM equation in its general form: 



for x G M, t > 0, where a, (3, 7, 8 are positive real numbers, [4]. The finite volume method 
is well known for its accuracy, efficiency, robustness and excellent local conservative prop- 
erties. Most often this method is employed to approximate solutions to hyperbolic con- 
servation laws. The system of Nonlinear Shallow Water Equations (NSWE) is a classical 
example of the successful application of modern finite volume schemes to water wave prob- 
lems. 

A wide range of numerical methods have been employed to compute approximate solu- 
tions to dispersive wave equations of KdV-BBM type : finite difference schemes [10, 50], 
finite element methods [9, 33, 2] and spectral methods [38, 39, 15, 35]. Recently discontin- 
uous Galerkin schemes have also been employed to dispersive wave equations [51, 30, 18], 
(the list is far from being exhaustive). However, the application of finite volume or hybrid 
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FV/FD methods remain most infrequent for this type of problems. To our knowledge, only 
a few recent works are in this direction [3, 17, 6, 47, 42, 12]. 

In order to apply the finite volume method to the KdV-BBM equation (1.1), we rewrite 
it in a conservative form, including a nontrivial evolution operator, an advective and a 
dispersive flux functions. In the finite volume literature there exist several ways to approx- 
imate these fluxes. For the advective part we test three different numerical fluxes each one 
representing a particular family of finite volume method: 

• average flux (m-scheme), 

• central flux, (KT-scheme) as a representative of central schemes, [34, 29], 

• characteristic flux (CF-scheme), as a representative of upwind schemes and lin- 
earized Riemann solvers, [20, 22]. 

The dispersive term is discretized using simply the average flux, while high order approx- 
imations are used for the BBM term (ju xx t). The central flux and the characteristic flux 
are widely used in the case of conservation laws. On the other hand the average flux, 
known to be unstable for conservation laws, performs equally well. 

The evaluation of the numerical flux functions require approximate values of the solution 
at the cell interfaces. The order of the approximation determines the space accuracy of the 
underlying finite volume scheme. We consider first order, taking simply piecewise constant 
approximations, as well as high order schemes. The high order accuracy is achieved through 
application of various reconstruction techniques such as TVD [46], UNO [26] and WENO 
[31]. 

The time discretization of (1.1) is based on Runge-Kutta methods. The stability of the 
resulting system of ode's depends on the interplay between the BBM term (7 u xxt ) and the 
KdV type dispersive term (6 u xxx ). An explicit discretization of the ode system is sufficient 
when these terms are of the same order. Thus, Strong Stability Preserving Runge-Kutta 
(SSP-RK) methods, which preserve the TVD property of the finite volume scheme, [44, 24] 
are used for the explicit discretization. 

However, when 7 C (5 the resulting semidiscrete system of ode's is highly stiff and 
therefore implicit methods with strong stability characteristics are preferable. To balance 
the high computational cost of fully implicit methods and stability considerations we rely 
on Implicit-Explicit Runge-Kutta (IMEX) methods, [1]. Indeed IMEX RK methods turned 
out to be well suited for the time discretization of the KdV-BBM equation (1.1) exhibiting 
excellent stability behavior. 

The validated numerical method is applied to study the KdV-BBM equation (1.1) in a 
systematic way through a series of numerical experiments. In particular, we focus on the 
following issues: 

• accuracy of the finite volume method for solitary wave propagation and invariants 
conservation 

• dispersive shock formation (we underline that the finite element as well as spectral 
methods break down for this experiment while the finite volume method provides 
robust and accurate results) 

• interactions of solitary waves (overtaking collisions) 
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The paper is organized as follows. In Section 2 the governing equation (1.1) is presented 
briefly along with its basic properties. In Section 3 the finite volume discretization as well 
as fully discrete schemes are presented in details. In Section 4 we validate the discretization 
procedure by comparisons with analytical solution. Several important test cases are also 
presented. 

2. Dispersive water wave model equation 

We present briefly the mathematical model under consideration and some of its basic 
properties. The KdV-BBM equation takes the following general form: 

ut + au x + (3uu x — ju xxt + 8u xxx — 0, (2.1) 

where x G M, t > 0, u denotes the free surface elevation above the still water level u = and 
a,/3,7, S are positive real numbers. Equation (2.1) incorporates nonlinear and dispersive 
effects and has been suggested as a model for surface water waves in a uniform channel 
with flat bottom, cf. ([4, 19]). 

When 5 = 0, (2.1) reduces to the BBM equation [4], while taking 7 = leads the 
celebrated KdV equation [28]. The KdV-BBM model (2.1) has been studied thoroughly in 
the past and the Cauchy problem is known to be well-posed in appropriate Sobolev spaces, 
at least locally in time. Also the well-posedness of some initial-boundary value problems, 
including the initial-periodic boundary value problem, can be proved, cf. e.g. [4, 7, 19] 
and the references therein. 

One may easily check that equation (2.1) admits exact solitary wave solutions of the 
form: 

u(x,t) = 3^psech 2 QJ^L^L ( x _ C A , (2.2) 

that travel rightwards with a given speed c s . We are going to exploit this solution below 
in order to validate our discretization procedure and measure the order of convergence of 
proposed numerical schemes. Further it is well known that (2.1) possesses two quantities 
invariant under its evolution dynamics. Assuming either the solution has compact support 
or u — y x — y ±00, one can easily check that quantities 

h(t) = / u(x,t)dx, h(t) = / (u 2 (x, t) +ju 2 x (x,t)) dx, (2.3) 
Jm. Jr 

are conserved in time, i.e. Ii{t) = /i(0), ^(t) = -^(0), V£ > 0. The invariant I\ reflects the 
physical property of the mass conservation, while invariant I2 can be assimilated to the 
generalized kinetic energy. Invariants conservation is a fundamental property important 
not only for theoretical investigations but also for numerics since it allows to validate 
numerical schemes and to quantify the accuracy of the obtained results. 

For more realistic situation one has to consider bidirectional models with uniform or 
variable bathymetry cf. e.g. [8, 40]. For a systematic numerical study of such Boussinesq 
type systems using finite volume methods analogous to those presented in this paper, 
including the runup algorithm we refer to [16]. 
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3. Finite volume discretization 



We proceed to the discretization of (2.1) by a finite volume method. Our motivation 
stems from the observation that the KdV-BBM equation can be seen as a dispersive per- 
turbation 1 of the following inviscid Burgers equation: 

u t + (au + t^u 2 ) x = 0. 

Consequently, the proposed finite volume schemes are based on the corresponding schemes 
for scalar conservation laws. A special treatment is introduced for the discretization of 
dispersive terms. 

Let T = {x{}, i G Z be a partition of R into cells Cj = (Xj_i , x i+ i), where Xi = 
+ x i _i)/2 denotes the midpoint of the cell Cj. Let Axi = x i+ i — x^i denote the 
length of the cell Cj and let Ax i+ i = x^+i — Xj. Herein, we assume the partition T to be 
uniform, i.e. Axi = Ax i+ i = Ax, i £ Z. For a scalar function w(x,t) let Wi denotes its 
cell average on Cf. 

1 f 

Wi(t) = — — / w(x,t)dx. 
Ax J c . 

We rewrite (2.1) in a conservative-like form: 

(/ - -yd&ut + [F(u)] x + [G(u xx )] x = 0, (3.1) 

where the advective flux is F{u) = au + | w 2 and the dispersive flux is G(v ) = 5v. We 
underline that F is a convex flux function. A simple integration of (3.1) over a cell Cj 
yields: 

d 



dt 



+ 



Ui(t) - (u x (x i+ i,t) -u^x^ i,t) J 
F(u(x i+ i,t))-F(u(Xi_i,t)) 



Ax L 



+ 



Ax 



G(u xx (x i+ i,t)) - G(u xx (xi_i,t)) 



(3.2) 

where the values of the advective and dispersive fluxes on the cell interfaces have to be 
properly defined. 

3.1. Semidiscrete scheme. We proceed to the construction of the semidiscrete finite 
volume approximation. Let xc t be the characteristic function of the cell Cj. We define a 
piecewise constant function Uh(x,t) = J2iez Ui(t)xcX x )i where Ui(t) are solutions of the 
following system of ordinary differential equations: 



d 

dt. 



Ui - 



7 

Ax 



U i+1 - 2Uj + Uj-i 
Ax 



Ax 



Ax 



^x =0, (3.3) 



with initial conditions defined as a projection onto the space of piecewise constant functions 
on T: 



^(0) = — / u(x, 0) dx, i G Z. 
Ax J a 



^Since the wave is assumed to be weakly nonlinear and weakly dispersive. 
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In (3.3) J 7 and Q denote the advective and the (KdV-type) dispersive numerical fluxes 
respectively. More specifically, J- i+ i = J~(U^i,U^_i) and G i+ i = Q(Wf^i,W^_i) are 

approximations of F(u(x i+ i,t)) and G(u xx (x i+ i , t)) respectively at cell interfaces. Values 
U^ + i, U^i are approximations to the point value u(x i+ i , t) from cells C iy C i+ i respectively, 
while W L i and W^i are corresponding approximations to the point value of the second 
derivative u xx (x i+ i,t). All quantities U^ +l , U^i as well as W^_i, Wf^i are computed by 
a reconstruction process described below (see Section 3.1.2). 

3.1.1. Advective and dispersive numerical fluxes. Over the last twenty years numerous nu- 
merical fluxes J 7 have been proposed to discretize advective operators [41, 25, 37, 23, 5]. 
We select three quite different flux functions. Namely, we consider a simple average flux 
J= m , a central type flux J= KJ \ [29, 34] and a characteristic flux T CF ',[20, 21, 22] : 

r»{U, V) = F (¥±X.) , (3.4) 



F KT (U, V)= 1 - {[F(U) + F(V)) - A(U, V) [V - U)} , (3.5) 
F CF (U, V) = \ {[F(U) + F(V)\ - A(U, V) [F(V) - F(U)]} . (3.6) 

The average flux is perhaps the simplest one and is known to be unconditionally unsta- 
ble for nonlinear conservation laws. However, this flux shows very good performance for 
dispersive waves (see Section 4). 

The central flux is of Lax-Friedrichs type and is a representative of the family of central 
schemes. The operator A in the KT-scheme is related to characteristic speeds of the flow 
and is given by this expression: 

A(U,V) = m a x[\F'(U)l\F'(V)\]. (3.7) 

The characteristic flux function is somehow similar to the Roe scheme [41] and the 
operator A in this case is defined as: 



A(U, V) = sign (f'(^±-^) J = sign (a + /? 



For the dispersive numerical flux Q we choose to work with the average flux function 
(3.4): 

W + R 

g(W,R) = 5^—, (3.9) 

where W and R are standard central approximations of the second derivative from each 
side. The numerical flux Q can be evaluated either using simple cell averages, denoted by 
Q m , or higher order approximation based on a reconstruction procedure, denoted by Q lm . 
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3.1.2. Reconstruction process. The values U^ +1 , U^_ L are approximations to u(x i+ i,t) from 

cells Cj and Cj+i respectively. The simplest choice is to take the piecewise constant ap- 
proximation in each cell: 

Uf + , = U u U* +1 _ = U i+l . (3.10) 

The resulting semidiscrete finite volume scheme is formally first order accurate in space. 
To achieve a higher order accuracy in space, we have to adopt more elaborated reconstruc- 
tion process. The main idea is to use the cell averages U to reconstruct more accurate 
approximation to the solution at cell interfaces u(x i+ i,t). For this purpose we consider 
three different reconstruction methods: the classical MUSCL type (TVD2) piecewise lin- 
ear reconstruction [27, 48], the UN02 reconstruction [26] and WENO type reconstructions, 
[31]. 

• The classical TVD2 scheme uses a linear reconstruction : 
U L + i =U + \cf>(n)(U i+1 - Ui), U*, = U i+1 - ±</>{r i+1 )(U i+2 - U i+1 ), (3.11) 

where = ^ l ~ u '~f , and 6 is an appropriate slope limiter function, [46]. There exist 
many possible choices of the slope limiter. Some of the usual choices are 

- MinMod (MM) limiter : (f)(9) = max(0, min(l, 9)), 

- VanLeer (VL) limiter : ^ - 



1+|0|' 

- Monotonized Central (MC) limiter : (f)(6) = max(0, min((l + 9)/2, 2, 29)), 

- Van Albada (VA) limiter : (f)(9) = f±g. 

The last three limiters have been shown to produce sharper resolution of discontinu- 
ities, and in our case less dissipative numerical results. The TVD2 reconstruction is 
formally second order accurate except at local extrema where it reduces to the first 
order. Reconstructions considered below were proposed to remove this shortcoming. 
The UN02, like the TVD2, is also a linear reconstruction process which is second 
order accurate even at local extrema. The values U , i , U , i are defined as 

U^ = Ui + ^Si, U^ = U i+1 -^Si +1 , (3.12) 

where 

Si = m(S+,S7), <?± = d i±k U T\D i±h U, 
d i+ iU = U i+1 - U u D i+ iU = m(DiU, D i+1 U), 

DiU = U i+l - 2Ui + m(x,y) = ^(sign(x) + sign(y)) min(|x|, \y\) 

The UN02 reconstruction is formally second accurate even at local extrema. 
We also consider WENO type reconstructions [31, 43]. Namely, we implement the 
3rd and 5th order accurate WENO methods, hereafter referred to as WEN03 and 
WEN05 respectively. For the sake of clarity, we present here only WEN03 scheme. 
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First of all we compute the 3rd order reconstructed values: 

uf\ = lm - u i+1 ), uf\ = hui-i + Ut). 

1 2 I 1 2 Z 

Then, we define the smoothness indicators: 

A) = (u i+1 - u t ) 2 , Pi = {Ui - u^) 2 , 

and constants do — |, d\ — g, do = di, d\ = do- The weights are defined as: 

a,Q «o olq dii 

COq = , CO 1 = , COq = — — , CO 1 = — — , 

«0 + Oil CiQ + Oti «o + Oil a + Oil 

where = <Sj = and e is a small, positive number (in our computations 
we set e = 10~ 15 ). 

Finally, the reconstructed values are given by formulas: 

Remark 1. JTie elliptic operator approximation in (3.3) is orziy second order accurate. In 
the case where a high order WENO reconstruction is used, we need to increase also the 
elliptic solver accuracy. For example, the following semidiscrete scheme: 



d_ 

It 



Ui-! + lWj + Uj+i Uj+i - 2Uj + Uj.i 
12 7 Ax 2 



+ n lzl +im + n 1± i =Q (3l4) 



where % = — + — Gi-i) is a fourth order approximation. Thus in 

the WEN03 case a global third order accuracy is observed, while for WEN05 interpolation, 
we profit only locally by the 5th order accuracy of the reconstruction, cf. Section 4-1- 

Remark 2. In computation of the dispersive flux we distinguish between the simple aver- 
aging of cell centered values in Q m and of Q lm , where higher order reconstructions of the 
second order derivatives are used. 

3.2. Fully discrete schemes. We consider now fully discrete schemes for the ode system 
(3.3). The time discretization is based on Runge-Kutta type methods. Explicit schemes 
based on TVD preserving RK-methods are presented. In certain cases where stiffness 
becomes dominant, we use an implicit-explicit strategy based on IMEX type RK-methods. 

3.2.1. Explicit schemes. The initial value problem (3.3) can be discretized by various meth- 
ods. When the parameter 7 is of the same order as 5 the system of ode's appeared to be 
non-stiff and therefore can be integrated numerically by any explicit time-stepping method. 
We use a special class of Runge-Kutta methods that preserve the TVD property of the 
finite volume scheme, [44, 24, 45]. 
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Let At be the temporal stepsize and let t n+l = t n + At, n > be discrete time levels, 
then (3.3) is an initial value problem of the form 



TU' = L(U), 



(3.15) 



where U = {Ui}, i G Z, T = I + [—7, 27, —7] /Ax 2 is a tridiagonal matrix and L is a 
nonlinear operator incorporating the contribution of the numerical fluxes J 7 , Q. Assuming 
at time t n , U n is known then U n+1 is defined by 



At 
Ax 



U n+1 = U" - ^ ^& i T _1 .L(U n ' i ), 



U nj = U n 



At 
Ax 



-1 



(3.16) 



where the set of constants A = (a^), & = (bi,...,b s ) define a s— stage Runge-Kutta 
method. The following tableau are examples of explicit TVD RK-methods which are of 
2nd and 3rd order respectively 












1 





1 


1 


1 




2 


2 

















1 








1 


1 

4 


1 

4 





1 

2 


1 


1 


2 




6 


6 


3 





(3.17) 



In our computations we mainly use the 3-stage third order method. 



3.2.2. Implicit- Explicit schemes. As the parameter 7 decreases to zero the semidiscretiza- 
tion of the KdV-BBM equation leads to a stiff system of ode's. To solve efficiently this 
system we apply an IMEX type RK- method, [1]. The linear dispersive terms are treated 
in an implicit way while the rest of the terms are treated explicitly. Numerical evidence 
shows that IMEX methods exhibit excellent stability and handle stiffness in an efficient 
and robust way even in the limiting case 7 = 0. 

We consider an s-stage Diagonally Implicit Runge-Kutta (DIRK) method, properly cho- 
sen, that is given by the tableau 

a u • • • T\ 

; (s.is) 

Qsl a s2 ' ' ' O-ss Ts 

h b 2 ■■■ b s 



A 


T 


b 
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and an s + 1 explicit Runge-Kutta method 



A 


T 


b 








• • 


■ 








a u 


•■ 


• 





fl 


021 


0,22 ' ' 


• 





T~2 


o s i 


a S 2 ■ ■ 


o ss 







&i 


b 2 ■■ 


• b s 








We rewrite system (3.15) in the form 

TU' = .F(U) + DU, 



(3.19) 



(3.20) 



where D is the five-diagonal matrix 5[—l/2, 1, 0, —1, l/2]/Aa; 3 coming from the discretiza- 
tion of the KdV term when we use the numerical flux function Q m . Then the fully discrete 
scheme can be written in the form 



-l 

+ AtaiiD)U® = TU n - At ^ a^J"(U^) - At ^ a, 3 DU (j) , 



3=1 
s 



TU 



71+1 



TU n - AtJ2 ^-^(U (i) ) - At ^ 6 i DU^' ) . 

3=1 3=1 



(3.21) 
(3.22) 



We employ four IMEX RK-methods of different number of stages, orders of accuracy and 
stability properties. In particular we consider the following pairs, [1] 

• A two stage third order DIRK method and a corresponding three stage, third order 
accurate ERK method with 7 = (3 + v / 3)/6. The resulting IMEX method is third 
order accurate. 



(3.23) 



7 





7 


1 - 2 7 


7 


1-7 


1 


1 




2 


2 

















7 








7 


1-7 


2(1-7) 





1-7 





i 

2 


1 

2 





A two stage second order DIRK method which is stiffly accurate, with 7 = (2 — 
V / 2)/2. The corresponding ERK is a three stage second order accurate method 
with § = —2-^/2/3. The resulting IMEX combination is second order accurate. 



(3.24) 





7 





7 


1 


-7 


7 


1 


1 


-7 


7 



















7 










7 


5 


1 - 


5 





1 





1 - 


7 


7 





A three stage third order DIRK stiffly accurate method with larger dissipative 
region than (3.24). The corresponding ERK is a three stage third order method. 
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Method 


At/ Ax < 


(3.23) 


1/4 


(3.24) 


1/5 


(3.25) 


1 


(3.26) 


1 



Table 1. Stability of IMEX for the KdV equation (a = (3 = 5 = 1, = 0) 



The resulting IMEX pair is third order accurate. 



0.4358665215 
0.2820667392 
1.208496649 





0.4358665215 
-0.644363171 






0.4358665215 



0.4358665215 
0.7179332608 
1 



1.208496649 -0.644363171 0.4358665215 


0.4358665215 

0.3212788860 0.3966543747 

-0.105858296 0.5529291479 0.5529291479 







1.208496649 



-0.644363171 0.4358665215 





0.4358665215 
0.7179332608 
1 



(3.25) 



A four stage, L-stable DIRK method with rational coefficients. The corresponding 
ERK is a five stage third order method. The resulting IMEX method is third order. 
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2 
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3 
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7 
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1 


4 


4 







(3.26) 



We tested these IMEX methods in the case of the KdV equation with a = (3 = S = 1, 
7 = 0. In Table 1, we summarize the constraints for the timestep At, purely in term of Ax, 
to obtain a stable solution. IMEX methods (3.25) and (3.26) exhibit excellent stability 
behavior. 



4. Numerical results 

In this section we present a series of numerical results aiming to show the performance 
and robustness of discretization procedures described above. There are many possible 
combinations of numerical fluxes, types of reconstruction and slope limiter functions. We 
begin by examining the accuracy of the methods by measuring the convergence rates in 
Section 4.1 and the preservation of the invariants in Section 4.2. The ability of the schemes 
to capture a solitary wave solution is demonstrated in Section 4.3. Solitary wave collisions 
are studied in Section 4.4. Finally, a dispersive shock wave formation is investigated in 
Section 4.5. 
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(a) UN02 MinMod 



(b) WEN03 



Ax 


Rate(El) 


Rate(££°) 




Ax 


Rate(El) 


Rate(££°) 


0.5 


2.000 


2.015 




0.5 


2.604 


2.561 


0.25 


2.001 


2.014 




0.25 


2.790 


2.810 


0.125 


2.001 


2.012 




0.125 


2.905 


2.913 


0.0625 


2.001 


2.010 




0.0625 


2.974 


2.981 


0.03125 


2.001 


2.008 




0.03125 


2.968 


2.995 



Table 2. Rates of convergence : CF-flux 



Remark 3. The solution of the linear system involved in (3.15) and (3.21) is obtained 
by a variation of Gauss elimination for tridiagonal systems with computational complexity 
O(d), d— being the dimension of the system. 

4.1. Rates of convergences, accuracy test. We consider an initial value problem for 
(2.1) with periodic boundary conditions in [—100,100]. We take for simplicity a — /3 = 
7 = 5 = 1 and consider a solitary wave solution of the form (2.2) with c s = 1.1. We take a 
uniform mesh h = Ax = 200 /N and compute the solution up to T = 100 using the three 
stage third order explicit SSP-RK method (3.17) with time step At = T/M. The errors 
are measured using the discrete scaled norms E\ and E%°, [30] 



E 2 h (k) - 



\U k \ 



\U k \ 



\u°\ 



\u k \ 




h,oo/ 



\u°\ 



h,oo: 



\u k \ 



h,oo 



max I U- 



i=l,...JV 



where U h = {U k }f =1 denotes the solution of the fully-discrete scheme (3.16) at the time 
t k = k At. The numerical rate of convergence is defined by 

Rate = l °s(E hl /E h2 ) ^ 
log (h 1 /h 2 ) 

for two different mesh sizes hi,h 2 . 

We perform several tests using the TVD2, UN02 and WEN03 reconstructions. Nu- 
merical solutions are computed with CF, KT or average fluxes. Table 2 shows the rates 
of convergence for the CF-scheme along with UN02 and WEN03 reconstructions. We 
observe the theoretical 2nd order convergence for the average, TVD2 (not reported) and 
UN02 schemes. The WEN03 reconstruction in conjunction with improved elliptic inver- 
sion scheme (3.14) gives us the expected 3rd order convergence. Rates in Table 2 are 
obtained with the most dissipative MinMod limiter function, while other limiters yield 
slightly sharper results. Moreover, the convergence results for the average m— flux and the 
KT numerical flux are qualitatively identical to those of CF. Analogous convergence rates 
were obtained using the IMEX methods. 
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4.2. Invariants preservation. As already mentioned in Section 2, (2.1) admits at least 
two quantities (2.3) which remain constant under the equation dynamics. We investigate 
the conservation of these quantities by computing their discrete counterparts: 



The observation of invariants during numerical computations (3.16) may also give an idea 
on the overall discretization accuracy. 

The initial value problem for (2.1) with periodic boundary conditions is considered. We 
seta = /3 = 7 = <5 = l and consider a solitary wave solution with celerity c s = 1.5. We 
compute its evolution up to T = 200 using Ax = 0.1 and At = Ax/2. 

The first observation is that the mass of the solitary wave l\ = 13.41640786499 is 
preserved in all computations independently from the choice either of the numerical flux, 
reconstruction method, or the slope limiter function. 

The behavior of 1% is quite different. Figure 1 shows the evolution of the solitary wave 
amplitude and of the invariant 1%. The numerical solution is obtained using J 7 " 1 , T CF and 
J-"^ T numerical fluxes along with TVD2 and UN02 reconstructions. The limiter MinMod 
is used and the dispersive flux is computed with Q lm flux function. The behavior of CF 
and KT schemes is almost identical. Perhaps, the CF-scheme is slightly less dissipative 
than the KT-scheme. However, the m-scheme appears to be the least dissipative. 

For both KT and CF fluxes, the TVD2 reconstruction preserves neither the invariant 
I2 nor the amplitude of the solitary wave. In the same time UN02 reconstruction shows 
excellent behavior. Despite its simplicity, the m-scheme, using J 7 " 1 and Q m , performs very 
well too in preserving and the solitary wave amplitude. 

In Figure 2 we show the influence of the dispersive flux Q m , Q lm choice. One observes 
that G lm flux shows better behavior than the simpler G m flux. A comparable performance 
is achieved with CF-scheme using WEN03 and WEN05 reconstructions. 

Finally, in Figure 3 we show a comparison between the various slope limiter functions 
(Minmod, Van Albada, Van Leer and MC) tested with CF-scheme. MinMod limiter ex- 
hibits a small dissipative effect, while other limiters we tested show comparable behavior. 
The choice of the time-stepping method do not induce any difference. 

4.3. Propagation of solitary waves. We continue the presentation of numerical results 
by the classical test-case of a solitary wave propagation. This class of solutions (2.2) plays 
a very important role in the nonlinear physics and any practical numerical scheme should 
be able to compute with good accuracy this type of solutions. For simplicity, we will set 
to unity all coefficients a = /3 = 7 = 5 = 1 in (2.1). 

A large- amplitude solitary wave travels rightwards with the speed c s = 1.5. Its propa- 
gation is computed up to T = 100 with discretization parameters Ax = At = 0.1 using 
KT and CF numerical fluxes and TVD2 reconstruction. In both cases we use the Van 
Albada limiter. In Figure 4 we compare the analytical solution with the numerical one. 
Figure 4(b) is a magnification of the solitary pulse showing that the solitary wave shape is 
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ference on the vertical axis with respect to Figure 1). 



perfectly retained. Also we note that up to the graphical resolution, all curves are undis- 
tinguishable. In order to observe the differences between these solutions we present in 4(c) 
the error Eg = log 10 \u exact (x, 100) — U(x, 100) |. This shows that the difference between the 
numerical and the exact solution is analogous in all the cases and very small. 

The behavior of the numerical solutions can be better understood by analyzing the so- 
called effective equation, that is the p.d.e that the numerical scheme satisfies up to the 
order of the method. Obtaining an effective equation is not always feasible. In the case 
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construction : 'V': Minmod, '0': MC, Van Albada, V: Van Leer. 
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of the m-scheme for the KdV-BBM equation (2.1), the numerical solution Uh satisfies the 
following effective equation: 



Uh,t + ctu h .x + f3u h u htX - ju h)XXt + 5u Kxxx 

a 2 (- d H l £ -J. 

T L\X I Ufi,xxx ~r „ Mfttt/j T . ">h,x"'h,xx > j^h.xxxx ~ ^h,xxxxt 

\ o o 4 4 12 



0. (4.2) 



On Figure 5 we illustrate some artifacts of the numerical discretization for the pure BBM 
equation (5 = 0). In Figure 5(a) one can observe a small dispersive tail coming mainly 
from nonlinear terms discretization. The amplitude of the tail is related to the order of the 
method. Taking Ax ten times smaller leads the reduction of the amplitude by two orders 
of magnitude, as it can be observed on Figure 5(b). The explanation of these phenomena 
is contained in the straightforward analysis of the effective equation (4.2). 

We underline that the smallest tail is produced by the m-scheme and the largest by the 
KT-scheme. This shortcoming can be further reduced by UN02 or WEN03 reconstruction 
procedures. We conclude that a detailed study of solitary wave interactions would require 
a combination of a higher order method with a finer grid resolution. 

4.4. Solitary wave overtaking collisions. The solitary wave solutions (also known as 
solitons) of the celebrated KdV equation (a = /3 = <5 = l,7 = 0) have a well-known 
property to interact in an elastic way during an overtaking collision. In other words, the 
solitary waves retain their initial shape after the interaction, cf. [14]. Contrary to the KdV 
equation, the overtaking collision of two solitary waves of the BBM model and in general 
of the KdV-BBM equation is not elastic. Interacting solitary waves change in shape and 
also a small dispersive tail appears after the process. However, a nonlinear phase shift can 
be still observed even in the KdV-BBM equation. 
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Figure 4. Comparison between the analytical and numerical solutions: 
analytical solution, — : CF-TVD2, - -: KT-TVD2, -.-: m-scheme. 



Here we study the overtaking collision of two solitary waves of the KdV-BBM equation 
with a — /3 = 7 = 8 — 1. Solitary waves are located initially at X\ = —50 and Xi = 50 
with speeds c s = 1.5 and c s = 1.1 respectively. At t = we have two well separated pulses 
and the wave behind (left) propagates faster. Space and time variables are discretized 
with Ax = At = 0.01 to capture this process accurately. The solution is computed using 
the CF-scheme and three types of reconstruction: TVD2 with Van Albada limiter, UN02 
reconstruction with MinMod limiter and WEN03 method, and with the third order explicit 
SSP-RK method. 

The invariant I{* = 18.915498698 is conserved with the digits shown in all cases. With 
the invariant l\ the situation is slightly different: UN02 and WEN03 schemes preserved 
the value 1^ = 15.0633, while the more dissipative TVD2 reconstruction yields 1% = 15.063. 
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Figure 5. Dispersive artifacts of the equivalent equation: . . . : Analytical 
solution, — : FCVF-TVD2, - -:KT-TVD2, -.-:FCVF-WEN03 

Figure 6 shows the interaction process at several time instances in the left column, while 
the right column shows the corresponding magnification of the dispersive tail. Essentially 
no difference can be observed among various numerical solutions even in the magnified 
region, up to the graphical resolution. Additional snapshots aiming to illustrate the in- 
teraction process are shown on Figure 7. We observe that the solitary waves propagate 
connected as a single pulse with a single maximum for a small time interval contrary to 
bidirectional models [16] and to Euler equations (cf. [13]). 

Figure 8 shows the "elastic" collision of two solitons of the KdV equation (a — f3 — 5 — 1, 
7 = 0) up to t = 600. In this experiment we took Ax = At = 0.01 and 0.005 using IMEX 
method (3.25). Contrary to the analogous collision in the case of the BBM equation, we 
do not observe any new dispersive tails. Further magnification of the images show small 
artifacts of the order C(10" 6 ). The invariants are I\ = 12.280014566440 and If = 9.244 
for all the computations with Ax = 0.01. When a finer grid is considered, Ax = 0.005 
we do not observe any improvement in the conservation of the invariant I± while if was 
9.2442. Analogous conservation properties observed when we studied the collision for the 
KdV-BBM equation with the IMEX method (3.25) we observed that 1\ = 18.915498698945 
but no other improvement in the invariant if = 15.0633. 

4.5. Dispersive shock formation. It was proven that smooth solutions to the KdV 
equation tend to become highly oscillatory as the parameter 5 tends to zero, cf. [49]. 
These oscillatory solutions are sometimes referred to in the literature as dispersive shock 
waves. In this section we study numerically this special class of solutions. Recently, 
a discontinuous Galerkin method was employed to study the same problem [51] in the 
classical setting of the KdV equation. 

Namely we consider the KdV-BBM equation with a = — 1, 7 = 10 -5 and 5 = 0. A 
solitary wave solution (2.2) is taken as an initial condition with parameters a — /3 — 7 = 1, 
5 = and c, = 1.3. We underline that this initial condition is not an exact solution to 
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Figure 6. Inelastic overtaking collision of two solitary waves for the KdV- 
BBM equation 
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Figure 7. Inelastic overtaking collision of two solitary waves for the KdV- 
BBM equation (detailed view) 



the BBM equation under consideration, since the coefficient 7 is different. A fine grid 
with Ax = 0.001 is required to observe this phenomenon. We note that even much more 
accurate schemes [51] require almost the same resolution. Figure 9 shows the formation of 
a dispersive shock wave. The numerical solution is computed with four different methods: 
the m-scheme and CF-scheme with TVD2, UN02 and WEN05 reconstructions. The KT 
flux was also tested, producing almost identical to that of the CF-scheme. In all the cases 
we took At = Ax/10 except in the case of the WEN05 reconstruction where At = Ax/2. 

The invariant 7f = 7.493997530 conserving the digits shown during all simulations for 
all numerical schemes we tested. The behavior of 7^ is considerably different. Figure 10 
(left) shows that from the time the dispersive shock was formed, all numerical schemes, 
except the m-scheme, loose the conservation of the invariant 7^. As for the m-scheme the 
1 2 invariant was conserved to one decimal digit, during the whole simulation, see Figure 
10(b). 

On the other hand, when a solitary wave solution evolves for longer time intervals, using 
for example the m-scheme, we observe that solitary-wave-like structures are formed, cf. 
Figure 11, while retaining the conservation of the invariant 7^ up to one digit. Analogous 
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FIGURE 8. Elastic overtaking collision of two solitary waves computed with 
the KdV equation 
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Figure 10. Evolution of If 



behavior is observed for the KdV equation where general initial conditions evolved into 
series of solitary waves, cf. [14]. 

In Figure 12 we present the same experiment for the KdV equation (a = (3 = 1, 7 = 0, 
5 = 10~ 5 ) where the time integration is performed with the IMEX method (3.25) up to 
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Figure 12. Near the zero dispersion limit, KdV equation 
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T = 20 with discretization parameters Ax = 0.001 and At = Ax/2. We observe that 
the invariant 1^ is conserved with slightly, less accuracy while the 1^ = 6.572670686045. 
When we use the IMEX method (3.25) and the m-scheme in the case of the BBM equation 
we observe that the invariant Ig conserves 2 digits if = 4.49, while l\ = 7.493997530374 
conserving the digits shown. Thus we conclude that in this experiment (as also observed 
in all previous ones) the use of the IMEX method might improve the conservation of mass. 

5. Conclusions 

The main scope of the present article is to extend the framework of finite volume methods 
to scalar unidirectional dispersive models. We chose the celebrated BBM-KdV equation 
(2.1) as an important representative model arising in the water wave theory and having all 
main features of dispersive wave equations. 

The BBM-KdV equation can be also viewed as a dispersive perturbation of the inviscid 
Burgers equation. Consequently, our method relies on classical finite volume schemes which 
discretize the advection operator. Then, a special treatment was proposed for the KdV- 
dispersion term, while the BBM-dispersion required an elliptic operator inversion per each 
time step, hence, providing a physical regularization to numerical solutions. We propose 
and implement also several methods to obtain high order accurate schemes. 

The proposed discretization procedure is validated by comparisons with an analytical 
solitary wave solution. The order of convergence is measured as well as invariant preser- 
vation is studied extensively. The numerical method is applied to several important test 
cases such as a solitary wave propagation and a dispersive shock formation. We make also 
use of proposed higher order extensions to study the overtaking solitary waves collision for 
the KdV-BBM equation. 

The extension to more realistic bi-directional wave propagation models such as Boussi- 
nesq type equations [40, 36, 32, 16]. 
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